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SUMMARY 


The objective of this research is to formulate and 
evaluate a new numerical method for simultaneously and effi- 
ciently coupling an external subsonic potential flow and an 
interior viscous flow such that the two flows match at an 
interfacing boundary. Both a panel method and a simple 
point compressible vortex model are used for the outer 
potential field. The interior flow solvers which were used 
are the Navier-Stokes and Euler codes of T. J. Coakley 
(NASA Ames) and the Euler code of A. Verhoff (MCAIR) . In 
order to test compatibility/ the panel method is coupled to 
the less expensive Euler codes since the coupling procedure 
is identical to the Navier-Stokes code. 

The results of this study show significant efficiency 
improvements can be obtained over the uncoupled approach. 
Results also indicate the outer potential flow is best 
represented by the simple point compressible vortex model. 
The panel method couples smoothly to Coakley' s implicit 

code but is numerically incompatible as coupled with 
MCAIR' s explicit Euler code. 

An improved Navier-Stokes code is under initial 

development which extends MCAIR' s Euler code to include the 

necessary viscous terms. Results are shown for an infinite 

length channel with one wavy periodic wall with and without 
laminar separation. Further development is recommended for 
airfoil applications with the point compressible vortex 
used for the outer potential flow region. 
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SECTION I 
INTRODUCTION 


The primary objective of this effort is to develop an 
efficient coupled viscous/inviscid flow solver capable of 
accurately analyzing airfoil flows with laminar separation. 
A numerical method is formulated and evaluated which 
couples an external subsonic potential flow and an interior 
viscous flow such that the two flows match at an interfac- 
ing boundary located outside the region where viscous 
effects are important (see Figure 1). An existing Navier- 
Stokes Equation solver is used to compute the flow field 
for the inner viscous region including flow values on a 
paneled boundary. The outer region is governed by the 
LaPlace equation where all flow quantities can be deter- 
mined from the velocities on the surface of the paneled 
boundary. The matching problem is to solve for all flow 
quantities inside and on the panel boundary using the 
Navier-Stokes flow solver with outer boundary values 

determined by the coupled panel method. 

Objective 

• Predict Effect of Laminar Separation on Subcritical Performance 

• Efficient Time Dependent Solutions at Realistic Reynolds Number 

Approach 

1) Navier-Stokes Equation Solver 

2) Paneled Interface Boundary 



Figure 1. Numerical Approach for Airfoils With Laminar Separation 

This contract is part of a larger joint effort as 
shown in Figure 2 to analyze airfoils undergoing laminar 
separation. McDonnell Aircraft Company (MCAIR) had 
developed under this contract a coupled potential/Navier- 
Stokes method and evaluated it for accuracy and efficiency. 
R. J. Margason of NASA LaRC served as technical monitor for 
this effort as well as for the complementary research of 
Drs. K. N. Ghia, U. Ghia and G. A. Osswald of the Univer- 
sity of Cincinnati (Reference 1). 



NASA-Langley 

• Technical Guidance 

• Diagnostic Wind Tunnel Testing 


CRAD 


NASA 

Grant 


MCAIR 

Merge Panel Method With Navier-Stokes | 
Evaluate Accuracy and Efficiency 


Cooperative 

IRAD 


University of Cincinnati 

• Unsteady Navier-Stokes 

• Small Scale Separation Bubbles 

• Drs. Ghia/Osswald 


Figure 2. NASA-Industry-University Cooperative Research Effort 
Analytical Method for Laminar Separation 


The work of B. R. Gilmer and D. R. Bristow, MCAIR 
(Reference 2) serves as motivation to the MCAIR effort. A 
fully coupled integral boundary layer /panel method was 
developed to analyze airfoils at high angle of attack under- 
going massive turbulent separation. Three representative 
solution sets are shown in Figure 3 for differing Reynolds 
Number . Results agree very well with test data (Reference 
3). The advantage of using a Navier-Stokes equation solver 
in place of a boundary layer method is to increase applica- 
bility to more general flows. Airfoils undergoing laminar 
separation with possible reattachment are of particular 
interest to this study. This report describes our effort 
to develop such a coupled approach . 

Section II describes the formulation of the necessary 
panel relations between outer boundary velocities and panel 
velocities. Several approximations are described based on 
potential flow theory. The accuracy of the panel method is 
addressed in relation to applicable locations of the panel 
boundary. 



Stalled Airfoil Analvsis Proaram (ref. 2 ) Experiment (ref. 3) Reynolds number 
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The coupling procedure is discussed in Section III . 
The interior flow solvers used in this study are the Navier- 
Stokes and Euler codes of T. J. Coakley, NASA Ames 
(References 4-6) and the Euler code of A. Verhoff, MCAIR 
(Reference 7). The panel method is coupled only to the 
Euler codes in order to reduce developmental cost. Coupled 
results are shown using Coakley' s Euler code which indicate 
a significant increase in accuracy over the uncoupled code. 
A numerical incompatibility results between the panel 

method and the explicit Euler code of MCAIR as presently 
coupled. An alternative to the present coupling procedure 
is discussed. 

A simplified outer potential field concept is pre- 
sented in Section IV. A simple point compressible vortex 
is used to model the effect of airfoil circulation on the 
outer boundary as suggested by J. L. Thomas, NASA LaRC 

(Reference 8). The strength of the vortex is based on the 
airfoil pressure integrated lift coefficient, thus being 

applicable to viscous as well as inviscid flows. Results 
are shown which indicate the current coupling procedure 
with the panel method shows no benefit over the simple 

point vortex model of an outer potential field. 
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The approach that appears most promising for the analy- 
sis of airfoils with laminar separation is to add viscous 
terms to the MCAIR explicit Euler method while using the 
point compressible vortex for the outer flow region. 
Preliminary development of a laminar Navier-Stokes code is 
presented in Section V. The formulation is based on 
Riemann Variables in local streamwise coordinates. The 
advantages of this approach are discussed and preliminary 
results are given for infinite length wavy wall channel 
flows with and without laminar separation. Further 
development is recommended for airfoil applications. 

All calculations for this study were performed on the 
NASA LaRC computing complex. Most flow calculations were 
performed on the CYBER 203 vector machine. 
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SECTION II 
THE PANEL METHOD 


In this section, the fundamentals of the panel method 
are explained, the accuracy of the panel method in reproduc- 
ing outer velocities computed by the Navier-Stokes flow 
solver is estimated, and the extent to which the outer 
potential region may approach the airfoil surface is deter- 
mined approximately. The matching process between the 
panel method and the interior flow solver will be described 
in the next section. 


The panels are interconnected along a chosen C-grid 
line and closed along a vertical grid line downstream as 
shown in Figure 4. The objective is to establish panel 
singularity strengths that nearly reproduce the velocity 
vectors on the panels. It is not generally possible to 
establish a potential field that will reproduce these veloc- 
ities since the panel boundary is influenced by viscous and 
rotational effects. Either set of tangential or normal vel- 
ocity components alone will define unique potential fields. 


The outer potential field may be defined by one of the 
three choices shown in Figure 4. The panel boundary value 
problem is then easily set up and solved using Morino type 
boundary conditions (Reference 9). Given N panels, it is 
sufficient to solve N-l control point equations along with 
an auxiliary relationship in order to define the potential 
outer region. 


4>i+i - <J>i “ o 


2 [(A 

j - 


a i+l. 
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An. . ) a; + (A y . , 
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i = 1,2, . . .N-l 


The Nth equation satisfies zero total mass flux through the 
panel boundary when the outer potential is defined by 
tangential panel velocity components. 


0 = E 


j 4s j°j 


Panel length is represented by A S. When the potential is 
defined by normal panel velocities, the Nth equation 
controls the circulation on the panel boundary. 
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{<£i n t} = {0} = [A<,]{a} + [A 7 ]{7} 

o { = Local Source Strength 
7 j = Local Vortex Strength 

Choice of Outer Potential Flowfields: 

1 . {-y} = { U,}; {a}= — [A a ] -1 * [A T ]{-y} (Defined by Tangential Velocities) 

2. {ff} = {U n }; { 7 } = - [A ^] -1 • [A ct ]{( 7 } (Defined by Normal Velocities) 

3. Some Combination of 1. and 2. 

Figure 4. Definition of the Outer Potential Flow Region by the Panel Method 


The influence coefficients are constructed by the 
methods given in Reference 10 and are of the form: 

Ut i = • [ D °ij + Dy ij 'ij 
u „i - * kij °j + Ey ij Y jJ 

where the "+" superscript represents potential region 
values exterior to the panel boundary. 

Substitution of the appropriate relations shown in Figure 4 
yields the combined influence coefficients (B^j, c ij)‘ 
Given the first choice of potential fields we have: 

“t! - : p <■>„ A 0 ' 1 R Y ). . + D Yij j Yj = I B ijYj 

Yj = E j Yj 

j 


- (E a A a L A Y ) 
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Any of the three choices of outer potential will 
exhibit a velocity mismatch between the two regions since 
the inner region is generally nonpotential . The potential 
region can only be expected to match one velocity component 
along the panel boundary. This problem can be eliminated 
if we can accept a non-potential field generated as 
follows . The outer tangential velocities are computed 
using tangential panel velocities to define the singulari- 
ties, whereas the outer normal velocities are computed 
using normal panel velocities to define their own singulari- 
ties. The resulting outer region is generally non- 

potential since the two potential fields are selectively 
combined using appropriate velocity components from each 
potential field. The advantage to this approach is a 
smooth match between velocities along the panel boundary. 

The influence coefficients based on panel tangential 
velocities are computed using linear vortex and linear 
source panel distributions, whereas those based on normal 
velocities use linear vortex and constant source distribu- 
tions. The influence coefficients need to be computed only 
once before the Navier Stokes iteration process begins . 
After each iteration the outer velocities can be updated 
using the above relations. Region 1 is matched with Region 
3 (refer to Figure 1) when the panel velocities determine 
the outer- velocities which in turn satisfy the Navier- 
Stokes boundary value problem. 

A complete coupling of panel boundary conditions with 
the Navier-Stokes computer code was not attempted at this 
time. Instead it is instructive to consider a comparison 
between the velocities on a desired grid line from a compu- 
ted Navier-Stokes solution with those same velocities recom- 
puted by the previous panel relations along a chosen 
paneled boundary. This allows us to observe the compati- 
bility of numerical approximations between the Navier- 
Stokes code and the panel method . 

A Navier-Stokes solution is computed using Coakley's 
code (Reference 6) for a NACA 0012 airfoil at two degrees 
angle of attack and freestream Mach Number of 0.3 given a 
160 by 50 C-grid. The velocities are recorded along every 
grid line surrounding the airfoil. Approximately every 
third point in each grid line is compared with recomputed 
values using the panel relations based on an inner paneled 
grid line. For each case the inner paneled grid line is 
located three cells inside the grid line of interest. In 

Figure 5, the horizontal and vertical velocities (U x , Uy) 
are compared for the 40th grid line from the airfoil. This 
grid line is located approximately one chord length from 
the airfoil surface and extends ten chord lengths down- 
stream from the trailing edge. The plotted velocities are 
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perturbations from freestream in order to magnify the devia- 
tion between Navier-Stokes and potential velocities. The 
comparison at this station between outer velocities is good 
with two notable exceptions. First, a significant discrep- 
ancy occurs in the U x perturbations along the wake of the 
airfoil. Secondly, a residual vertical velocity error ic 
apparent along the outer boundary. The overall comparison 
indicates the panel approximation can serve as a good bound- 
ary value approximation for the interior flow solver. 



NACA 0012 

M oo = 0.3 a = 2° Re = 10® 


Calculated Velocities on Outer Boundary 

— Navier-Stokes solution (ref. 6) 

— — — Panel method approximation (ref. 10) 
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Figure 5. Accuracy of the Panel Method Approximation 
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A comparison of the panel approximation for the ten 
recorded grid lines is shown in Figure 6. An RMS error is 
computed relating Navier-Stokes velocities to velocities of 
the panel method for each grid line and plotted versus the 
height between grid lines/ h. The error is computed as f 


RMS ERROR 


z 

i s TOT 


k - «;» 2 + 

pan 


(U 


- U + ) 
n n 

pan 


h 

i 


where S refers to arc length along a grid line. For 
example/ the RMS error for the 40th grid line (refered to 
in Figure 5) is given at the height-to-chord ration of two. 
As we approach the airfoil/ the error grows rapidly due to 
the significant increase in diffused vorticity and relative 
magnitude of the perturbation velocities. The relation 

between h and the total number of grid points in Region 1 

is a measure of efficiency since the amount of computa- 

tional work is nearly proportional to the number of grid 
points. Work done by the panel method is negligible. The 

figure indicates a reduction of nearly 25 percent of the 
workload may be achieved with a minimal loss of accuracy by 
use of the coupled approach. This may be accomplished by 
eliminating about two thousand grid points in the outermost 
region. 



Figure 6. A Study on the Effect of Accuracy and Efficiency by Varying the Extent 

of the Outer Potential Flow Region 
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SECTION III 

COUPLING THE PANEL METHOD TO AN INTERIOR FLOW SOLVER 


A detailed description of the coupled approach between 
the panel method and two different interior flow solvers is 
given in this section. The selected interior flow solvers 
are both Euler codes in an effort to reduce development 
cost. The actual coupling process is identical between the 
panel method and either the Euler or Navier-Stokes version 
of the same code since the outer field is considered 
inviscid. Results are given and the various anomalies of 
the approach are discussed . 

The coupled iteration cycle is shown in Figure 7. The 
entire process starts with preliminary setup of grid geom- 
etry, initial flowfield and calculation of panel influence 
coefficients described in the previous section. The itera- 
tion cycle then begins with an unmodified interior flow 
solver iteration (e.g. Navier-Stokes) using frees tream 
boundary conditions on the outer grid lines . Panel singu- 
larities are then computed along a preselected panel 
boundary. This step is trivial since only one set of singu- 
larities corresponding to either tangential or normal per- 
turbation panel velocities are needed to define the outer 
potential field. The velocities on the outer boundary of 
the discretized flow solver are then computed using the pre- 
calculated influence coefficients. Constant enthalpy and 
entropy conditions are used to determine the remaining two 
boundary values on the outer grid line. The iteration 
cycle repeats with continually modified outer boundary 
values until convergence. 

The Euler problem would be overspecified if all four 
boundary values were specified as Dirichlet conditions. 
Only those boundary values that correspond to incoming 
characteristic waves are specified as shown in Figure 8 for 
Coakley's code. Namely, three flow values are specified at 
an inflow boundary and only pressure is specified at an 
outflow boundary. This is consistent with the one dimen- 
sional approximation of the wave equations as described by 
Coakley (Reference 4). The remaining boundary values which 
correspond to outgoing characteristic waves must satisfy 
the compatibility relations . 
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Panel Boundary 

czzZZZZZ za 


Outer Boundary of 
Discretized Flow Solver 


Coupled Iteration Cycle 



Choice of Outer Potential Flowfields 


1 . Panel Tangential Velocities 

2. Panel Normal Velocities 

3. Some Combination of 1 and 2 


Figure 7. The Coupled Navier-Stokes/Panel Method Approach 



Prescribed 

Outflow 

(P) 


Figure 8. Coupling the Panel Method With Coakley’s Implicit Euler/Navier-Stokes Code 
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The coupled method was applied to the NACA 0012 air- 
foil case at Mach 0.3 and 10 degrees angle of attack. 
Several calculations of lift coefficients are compared in 
Figure 9 varying the location of the outer boundary. The 
outer, boundary radius is reduced simply by truncating the 
outer grid lines from the original 120 by 30 C-grid. This 
enforces uniformity in the grids near the airfoil among the 
various cases. The outer velocities were calculated from 
both sets of panel velocity components as discussed in the 
previous section. The same cases were recomputed with the 
original Coakley code which uses freestream prescribed 
boundary values. 



Figure 9. Comparison of Lift Coefficient for Various Outer Boundary Locations 


Results are shown for four grid sizes ranging from a 
far field radius of 2.5 chord lengths to 20. Each C-grid 
extends 10 chord lengths downstream from the airfoil. Both 
methods approach nearly the same lift coefficient as the 
outer boundary radius approaches infinity. The limiting 
lift coefficients should not be expected to be exactly the 
same since the backplane is located 10 chord lenths down- 
stream from the airfoil in all cases. The coupled approach 
gives a more accurate lift coefficient for each case. The 
difference is magnified as the outer boundary radius is 
reduced. A major concern at this point is the variance in 
lift coefficient for the coupled method. 


12 




One of the more apparent reasons why the lift varies 
for the coupled approach is that circulation is not con- 
served in the Euler code. Figure 10 compares the circula- 
tions expressed as local lift coefficients for each C-grid 
line connected by the vertical grid line far downstream. 
The circulation should not vary for subcritical inviscid 
flows, but there exists a numerical error for this case. 
Since the circulation varies with grid line, it is apparent 
the choice in location of the panel boundary affects the 
circulation imposed on the outer boundary. The calculated 
lift coefficients for varying panel location are compared 
in Figure 11, indicating a considerable discrepancy. 



Airfoil J (Grid Line Index) Panel Outer 

Boundary Boundary 


Figure 10. Circulation Variation on C-Grid Lines 
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1.3 


NACA 0012 Euler Solution 
Mach 0.3 a = 10° 

O 120x25 grid 



09 l I I I I I I 

' 0 1 2 3 4 5 6 

Number of Cells Between Boundaries 

Figure 11. The Effect of Panel Location on a Fixed Outer Boundary 


Since the outer potential region starts at the panel 
boundary, it is desirable to locate the panels near to the 
outer boundary. However, when the panels are one half cell 
width inside the outer boundary, a conflict between the 
panel relations and the compatibility relations is intro- 
duced. For example, on the downstream boundary the compati- 
bility relation prescribes a zeroeth order extrapolation to 
determine tangential velocity, whereas the panel relations 
may call for a significant change. This compatibility rela- 
tion apparently becomes less of a conflict as the panels 
are moved nearer to the airfoil. 

The above calculations were made using panel relations 
between tangential velocities on the two boundaries and sim- 
ilar panel relations between normal velocities. The result- 
ing outer flow is generally non-potential since a different 
potential field was used to compute each velocity compo- 
nent. Panel relationships based only on tangential panel 
velocities have been used, but results are similar to those 
shown earlier . 
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The panel method is similarly coupled to MCAIR's 
explicit Euler code (Reference 7) for a second comparison 
of the generalized coupled approach. The boundary condi- 
tions are easily implemented in this code because the char- 
acteristics approach with local streamwise coordinates 
clearly leads to a straightforward outer boundary condi- 
tion. The prescribed Dirichlet boundary values are shown 
in Figure 12 for this extended Riemann flow solver. The 
coupling procedure between panel method and Euler method is 
similar to that with the previous Euler code. 


Prescribed Inflow 



Prescribed 
Outflow 
(R = q - as) 


Figure 12. Coupling the Panel Method With Explicit Euler Code 


The coupling of the panel method to the explicit code, 
however, results in a neutral instability which apparently 
is controlled only by reducing the frequency of updated 
outer boundary values. Figure 13 shows a NACA 0012 solu- 
tion history in which the panel method was used every iter- 
ation to update the appropriate boundary values. This 
solution is clearly unacceptable since the circulation does 
not approach a steady state. The same features are present 
when the outer boundary values are updated every two 
hundred iterations, but the magnitude of the oscillations 
has decreased. Convergence to steady state is achieved 
only after five hundred or more iterations are performed 
before each update of the boundary values. The rate of con- 
vergence, however, is slowed after each implementation of 
updated boundary values. A comparison of the lift coeffi- 
cients between the coupled method and the original code 
were not obtained because of the unwieldy convergence of 
the coupled method. Other attempts to alleviate the prob- 
lem were unsuccessful including under -relaxation, panel 
boundary movement, and different methods of generating the 
outer potential field. The reason for this instability may 
be the incompatibility between the panel method where infor- 
mation travels instantaneously to every other point in the 
potential field and the explicit Euler code where informa- 
tion can only propagate one cell length per iteration. 
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Outer Boundary 


■Paneled Boundary 



NACA 0012 
Mach 0.63 a =2° 




Panel Application Every 200 Iterations 


Figure 13. Panel Method Coupled With Explicit Euler Code 

Demonstration of Convergence Instability 


One of the difficulties with the general coupling pro- 
cedure used in this study is associated with the multiple 
outer boundary concept. Another approach would be to form 
one outer boundary by locating the panel boundary along the 
outer grid boundary. The problem would be altered to deter- 
mine the appropriate prescribed values which best approxi- 
mate an outer potential field in a least squares sense. 
The outflow prescribed value, R, shown in Figure 12 can be 
determined simply from the local values of Q and S using a 
constant freestream total enthalpy. At inflow, Q may be 
calculated similarly from R and S. Entropy prescribed at 
inflow would simply be the freestream value which leaves 
the set of incoming flow angles, 0, as the only unknowns. 
Unfortunately, this approach requires the solution of a 
large full matrix problem every iteration unless this 
matrix is grossly approximated by its diagonal terms only. 
This particular approach is not recommended for efficiency 
reasons? however, an investigation of this approach could 
lead to a better understanding of a properly coupled 
method . 
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SECTION IV 

A SIMPLIFIED OUTER POTENTIAL REGION 


A viable alternative to coupling a panel method with a 
Navier-Stokes solver is to simplify the outer potential 
flow method. Since the lift coefficient varies consider- 

ably for the coupled method shown in Figure 9, it is desir- 
able to maintain a large outer boundary radius to obtain 
the most accurate results. At large radius we can simplify 
the outer potential field by modeling the effect of the air- 
foil on the outer boundary with a simple point compressible 
vortex (Reference 8) located at quarter chord position. 
This method is coupled to Coakley's airfoil code as shown 
in Figure 14. The strength of the vortex is calculated 
from freestream velocity and pressure integrated lift coef- 
ficient. The appropriate outer boundary values are calcu- 
lated from superposition of the point vortex in the 
freestream flow with auxiliary relations of freestream 
total enthalpy and entropy. 



Prescribed 

Outflow 

(P) 


Figure 14. Simplified Outer Potential Field 

Point Compressible Vortex 


The NACA 0012 case at Mach 0.3 and ten degrees angle 
of attack is revisited using the coupled point vortex 
approach. The same four outer boundary radius solutions 
are calculated and compared with earlier results in Figure 
15. The point vortex solutions are slightly favorable to 
the coupled panel method and significantly better than the 
original code results. The point compressible vortex is 
simple to incorporate in the Euler /Navier-Stokes code and 
requires little computational overhead . The coupled panel 
method requires storage of two influence coefficient 

matrices (N by N) and uses approximately two percent more 
CPU time per iteration for a given mesh size as compared to 
the original code. The coupled point vortex approach has 
produced better results than the coupled panel methods 
discussed previously. It is simple, reliable, and more 
accurate for near field outer boundaries. 
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R 

Figure 15. Comparison of Outer Potential Flow Region Methods 


The point compressible vortex approach works equally 
well with Coakley's Navier- Stokes solver. Several airfoil 
cases have successfully been computed with this approach. 
A transonic RAE 2822 airfoil case is shown in Figure 16 for 
comparison between the original boundary conditions and 
modified values using the present approach. The lift as 
expected is significantly higher indicating the same prob- 
able trend as shown earlier. The point compressible vortex 
approach appears most favorable for analyzing airfoil flows 
with laminar separation and should work well with any basic 
flow solver. 
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Figure 16. RAE 2822 Airfoil Navier-Stokes Solution 

Q-Omega Turbulence Model 160 x 50 Grid 
Mach 0.73 a = 2.78° 
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SECTION V 

IMPROVED NAVIER- STOKES DEVELOPMENT 


A laminar Navier- Stokes equation solver is under devel- 
opment in an effort to improve the accuracy and reliability 
in comparison to available codes . An approach that appears 
promising is to add viscous terms to the MCAIR explicit 
Euler method (Reference 7) currently being developed under 
MCAIR Independent Research and Development (IRAD). The key 
advantage to this approach is the simplicity of the diagon- 
alized equations, all in wave equation form. This approach 
has been extended to laminar Navier-Stokes flows and devel- 
opment has progressed to the point of analyzing infinite 
length channel flows. In this section the basic formula- 

tion and advantages of the approach will be described, and 
preliminary results will be shown. 

The basic formulation is composed of Riemann Variables 
(Q,R) in local streamwise coordinates (s,n) shown in Figure 
17. Each equation represents a wave form indicating the 
proper physical direction in which the four unknowns 
(Q,R,0,S) are propagating. The source terms (Z^) are 
easily derived and are generally small. The Riemann Equa- 
tions are derived by combining the streamwise momentum equa- 
tion and continuity equation (in streamwise coordinates) 
while making use of the state equation and ideal gas assump- 
tion. The normal momentum equation and entropy transport 
equation are already in wave form. 

3f, 3f, 

+ X:— — — Z: — Wave Equation 

dt 1 os 



y 


q = Velocity Magnitude 
s a = Speed of Sound 

5 = Entropy 

6 = Flow Angle 


Figure 17. Riemann Invariant Form of Navier-Stokes Equations in Local Streamwise Coordinates 
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There are several advantages to this form of solution 
approach which deserve attention. The equations are mini- 
mally coupled since they are in diagonal form. They are 
only coupled through the source terms of the wave equations 
and the eigenvalues. The characteristics are unambiguously 
defined using the local streamline coordinate system. For 
each equation there is one and only one eigenvalue repre- 
senting a single wave propagation process. A simple and 
efficient algorithm gives a time dependent relaxation to 
steady state for an arbitrary grid. The algorithm is fully 
vector izable with vector lengths nearly equal to the number 
of grid points, making extensive use of CYBER 200 series 
technology. The formulation has proved versatile (Ref- 
erence 7) with a host of one, two and three dimensional 
Euler solutions. Several solutions ranging from incompress- 
ible Mach Numbers to high supersonic speeds have been com- 
puted. Viscous terms are easily incorporated into the 
equation source terms. 

Development of a laminar Navier -Stokes solver was ini- 
tiated for the simple case of two parallel plates in motion 
with constant distance between them. The plates initially 
move in the same direction at the same speed, then one 
plate is slowed to rest. The problem is to find a steady 
state solution and compare to the analytical solution. A 
two dimensional equation solver was developed first to 
solve this one dimensional problem, then extended to chan- 
nel flows . 

The parallel plate problem was computed simply with 
ten computational cells of equal size between the plates. 
Boundary conditions are composed of the no slip condition 
and specified plate temperatures. A Reynolds Number of ten 
based on distance between plates and a Prandtl Number of 
0.7 is specified with an upper plate speed at 0.5 Mach Num- 
ber . Results were obtained in a few hundred iterations 

which nearly reproduce the analytical steady state solu- 
tion. Results using extended Riemann Variables for this 
case were first given in Reference 7. Although the code 
uses a nonconservative formulation , total mass is conserved 
to 0.01 percent accuracy for this crude discretization. 

A channel flow problem is created simply by perturbing 
the shape of the lower surface. An 11 by 11 grid discret- 
ization is shown in Figure 18, where the lower surface rep- 
resents one period of a sine wave with a ten percent 
channel restriction. The flow is periodic with similar 
boundary conditions as the parallel plate case. The ini- 
tial flow guess is a linear variation in velocity compo- 
nents from top to bottom whereas entropy and speed of sound 
are assumed constant. The streamline contours shown are 
generated after a few hundred iterations representing a 
simple attached flow solution. 
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Re = 10 
Pr = 0.7 
11x11 Grid 

Periodic Boundary Conditions 


10% Channel Restriction 



T= 1.0 
M = 0.50 



Streamline Contours 


7 

T= 1.0 

A A 


Figure 18. Navier-Stokes Test Problem - Infinite Length Channel With Attached Flow 


A separated flow case was calculated by doubling the 
amplitude of the lower surface sine wave and raising the 
Reynolds Number to 1000. A refined 41 by 41 grid is used 
to discretize the flow as shown in Figure 19. The lower 
surface streamlines are shown at the four time intervals 
given. The flow was initialized in the same manner as the 
attached flow case. Each time unit corresponds to the time 
required for the upper plate to move one half the wave- 
length of the lower sine wave. After one time unit the 

flow is clearly separated and reaches its steady state loca- 
tion after two time units . The separation bubble grows in 
size and strength to that shown after eight time units. 
This solution represents sixteen thousand constant time 
steps for each grid point. This developmental code is 

formally first order accurate in both space and time and 

requires about one CPU minute per 1000 iterations for this 

case on a CYBER 203 using full precision operations exclu- 
sively. 
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Lower Surface Streamlines 
— 0.005 < ^ < 0.005 


Figure 19. Navier-Stokes Test Problem - Infinite Length Channel With Laminar Flow Separation 


The primary goal of this development is to improve the 
accuracy and reliability over existing methods. The use of 
Riemann. Variables in local streamwise coordinates forms the 
mathematical system of equations which dictates the proper 
numerical modeling of the physical wave propagation 
process. The strong relationship between mathematical, 

numerical and physical modeling in this formulation is 
expected to supply the means for an accurate and reliable 
code. Proper boundary conditions are easily implemented as 
a result of this formulation. 

A secondary goal is to develop a very efficient code. 
The developmental code has recently improved in efficiency 
by a factor of about six over the timings given above. The 
code now has a better streamlined logic and uses half pre- 
cision (32 bit) operations throughout. Local variable time 
steps were also implemented which requires fewer iterations 
to achieve steady state at the expense of time accurate 
modeling. A further efficiency increase of an order of mag- 
nitude may be achieved using higher order methods, multi- 
grid and/or local grid embedding concepts. 
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Future plans are to develop a Navier-Stokes airfoil 
code by extending MCAIR's Euler code using information 
learned from this developmental effort. Extended Riemann 
Variables will be used for a better shock modeling 
capability. Both algebraic and two-equation turbulence 
eddy-viscosity models will be implemented. The point 
compressible vortex approach will be used for an improved 
outer boundary condition until a better concept is devised. 
Subcritical attached flows will be treated initially, 
followed by flows with leading edge separation. Ultimately 
this development would be used to analyze wing/body 
configurations followed by complex aircraft configurations. 
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SECTION VI 
CONCLUSIONS 


The panel method couples smoothly with the implicit 
Coakley code resulting in a significant accuracy increase 
over the uncoupled code for any given outer boundary loca- 
tion. Three anomalies of the coupled approach are cited as 
follows : The accuracy of the solution may significantly be 

effected by changing the relative proximity between panel 
boundary and outer boundary. Secondly, the outer region is 
not potential since only some of the boundary values can be 
specified using potential relations. Finally, an insta- 
bility results when the panel method is coupled with the 
explicit MCAIR code. This instability decays when the 
outer specified boundary values are updated after every few 
hundred iterations. 

The results of this study show significant efficiency 
improvements can be obtained over the uncoupled approach. 
Results also indicate the outer potential flow is best 
represented by the simple point compressible vortex model. 

Preliminary development of a Navier-Stokes equation 
solver using Riemann variables has shown promising results. 
Infinite length channel flows are computed for attached and 
separated flow cases with relative ease. The formulation 
using local streamwise coordinates is simple, accurate and 
efficient and is recommended to be further developed for 
airfoil applications. 
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